% Reference: <<宽视场望远镜大气非等晕性的数值分析>> ---- 何超兰, Dotoral Dissertation
% Reference: <<An iterative algorithm for the position of time-varying phase screens in
% the atmospheric turbulence>>, Wenyi Shao and Hao Xian.
% this function is used for generating the params for splitted layers
% in multi-layer atmospheric turbulence model
% Written by Lu benchu, Taiyuan University of Technology
% email: benchuul@163.com
function [H_m, h_m] = mutilayer_generate(L_syst,layers,Cn2,loop_num)
%%% INPUT PARAMS:
% L_system: propagation length in turbulent area [m]
% layers: number of splitted phase screen layers
% Cn2(z): refractive structure constant, notice that 'z' with dimension [km]
% loop_num : 最大迭代循环计算的次数
%%% OUTPUT PARAMS:
% H_m: boundaries of splitted atmospheric layers [m]
% h_m：locations of splitted phase screens [m]

H_max = L_syst;

% 停止更新h_m的阈值
val_threshold = 1e-26;

% 赋予一些初始值
% Initializing the location of screen layers
h_m = zeros(1, layers); % 全局变量
h_m_record = zeros(1, layers);
h_interval = H_max/(layers+1);
for i = 1:layers
    h_m(i) = i*h_interval;
end

% splitted layers
H_m = zeros(1,(layers+1));
H_m(1) = 0; H_m(end) = H_max;

% 湍流层边界初始值
for i = 2:layers % computing H_m(2:end-1)
    H_m(i) = (h_m(i-1)+h_m(i))/2;
end

% z: altitude [km]
% 'z' represents 'h' in here
% fun_itgl = @(h) Cn2(h).*(h_i - h).^(2/3);

if ~strcmp(class(Cn2),'function_handle') % Cn2 is a constant
    % 直接均分大气路径为M层，屏位置位于路径均等分处
else % Cn2(z)
    for loop = 1:loop_num
        if min(h_m_record) == 1
            disp(abs(val_curr));disp(loop-1);
            break;
        end
        
        for i = 1:layers
            % 用梯度下降法迭代更新湍流屏位置h_m
            % References:《数值分析》-----p222,李庆扬，第五版
            
            % 当前循环，当前湍流层的最佳位置的算式值
            h_i = h_m(i);
            H_bot = H_m(i); H_top = H_m(i+1);
            fun_itgl = @(h) Cn2(h*10^(-3)).*((h_i - h).^(2)).^(1/3);
            val_left  = integral(fun_itgl,H_bot,h_i);
            val_right = integral(fun_itgl,h_i,H_top);
            val_curr = val_left - val_right;
            % 如果val_curr接近于0， 则考虑停止该次更新
            % abs(val_curr) > val_threshold
            if abs(val_curr) > val_threshold
                % computing the numerical differentiation
                val_diff = numerical_diff(Cn2,i);
                % update the location of h_m(i)
                h_m_pre = h_i - ((loop_num-loop)/loop_num)*(val_curr/val_diff);
                if i == 1
                    h_min = 0;
                    h_max = h_m(i+1);
                elseif i == layers
                    h_min = h_m(i-1);
                    h_max = H_max;
                else
                    h_min = h_m(i-1);
                    h_max = h_m(i+1);
                end
                h_m(i) = max(min(h_m_pre,h_max),h_min);
                h_m_record(i) = 0;
                % update the values of H_m
                for j = 2:layers % computing H_m(2:end-1)
                    H_m(j) = (h_m(j-1)+h_m(j))/2;
                end
                
            else
                h_m_record(i) = 1;
            end
            
        end
        % end of the loop_num
    end
    
% Cn2 为 与高度相关的函数 的情况
end

% Standarizing the dimensions: [km] to [m]
%H_m = H_m * 10^(3);
%h_m = h_m * 10^(3);

% numerical differentiation
% -----------------------------------------
    function val_diff = numerical_diff(Cn2,i)
        % 自变量为 h_i, 增量为dh, i 为顺序layer数:
        dh = 1e-5;
        % 前向差分
        h_mae = h_m(i)-dh;
        [H_mae, H_ato] = H_with_dh(i,h_mae);
        fun = @(h) Cn2(h*10^(-3)).*((h_mae - h).^(2)).^(1/3);
        val_forward = integral(fun,H_mae,h_mae)-integral(fun,h_mae,H_ato);
        % 后向差分
        h_ato = h_m(i)+dh;
        [H_mae, H_ato] = H_with_dh(i,h_ato);
        fun = @(h) Cn2(h*10^(-3)).*((h_ato - h).^(2)).^(1/3);
        val_backward = integral(fun,H_mae,h_ato)-integral(fun,h_ato,H_ato);
        % 数值微分
        val_diff = (val_backward-val_forward)/(2*dh);
        
        function [H_mae, H_ato] = H_with_dh(i,h_dh)
            switch i
                case 1
                    H_mae = 0;
                    H_ato = (h_dh+h_m(i+1))/2;
                case layers
                    H_mae = (h_m(i-1)+h_dh)/2;
                    H_ato = H_max;
                otherwise
                    H_mae = (h_m(i-1)+h_dh)/2;
                    H_ato = (h_dh+h_m(i+1))/2;
            end
            % end of the function
        end
        % end of the function
    end


% end of the function
end

